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The  focus  of  this  investigation  is  on  the  formulation  and  validation  of  a  methodology  for  the  estimation  of  a 
stochastic  linear  modal  model  of  a  structure  from  measurements  of  a  few  of  its  natural  frequencies  and  mode  shapes 
on  a  few  nominally  identical  samples  of  the  structure.  The  basis  for  the  modal  model  is  composed  of  the  modes  of  an 
approximate  representation  of  the  structure,  e.g.,  a  nonupdated  or  preliminary  finite  element  model.  Furthermore, 
the  variability  or  uncertainty  in  the  structure  is  assumed  to  originate  from  stiffness  properties  (e.g.,  Young’s 
modulus,  boundary  conditions,  attachment  conditions)  so  that  the  mass  matrix  of  the  uncertain  linear  modal  model  is 
identity  but  the  corresponding  stiffness  matrix  is  random.  The  nonparametric  stochastic  modeling  approach  is 
adopted  here  for  the  representation  of  this  latter  matrix;  thus,  the  quantities  to  be  estimated  are  the  mean  stiffness 
matrix  and  the  uncertainty  level.  This  effort  is  accomplished  using  the  maximum  likelihood  framework  using  both 
natural  frequencies  and  mode  shapes  data.  The  successful  application  of  this  approach  to  data  from  the  Air  Force 
Institute  of  Technology  joined  wing  is  demonstrated. 


Nomenclature 

en  =  overall  lit  of  modal  model 

K  =  modal  model  mean  stiffness  matrix 

Kxx  =  covariance  matrix  of  realized  values 
L  =  likelihood  function 

Sj  =  coordinates  of  measurement  point  i 

sL  =  coordinates  of  laser  collimator 
x{r)  =  realized  values  in  rth  experiment 
affj  =  measured  values  of  modal  coefficients 
S  =  uncertainty  parameter 

X  =  free  parameter  of  statistical  distribution  of  random 
matrices 

fix  =  mean  vector  of  realized  values 
<j>  I  =  finite  element  modes 

=  mode  shapes  projected  basis 

I.  Introduction 

HE  U.S.  Air  Force  is  responding  to  new  global  threats  by 
conceiving  aircraft  with  diverse  and  revolutionary  capabilities. 
These  capabilities  have  led  to  trends  in  smaller  fleet  sizes,  as  well  as 
more  diverse  flight  regimes  and  expanded  performance  limits.  Such  a 
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confluence  of  trends  renders  the  design,  analysis,  and  certification 
(DAC)  of  a  prototype  even  more  vital  to  ensuring  an  efficient  long- 
lasting  fleet.  An  integral  part  to  the  DAC  process  is  the  updating  of  a 
baseline  finite  element  model  (FEM)  with  ground  vibration  modal 
test  results  to  obtain  the  truth  model  for  a  prototype  structure.  This 
process  has  traditionally  been  a  deterministic  one,  yet  there  are 
indications  (e.g.,  [1])  of  a  definite  variability  in  the  structural 
properties  of  one  aircraft  component  vs  another. 

An  important  use  of  the  updated  FEM  has  been  the  construction  of 
a  reduced,  modal  model,  which  is  then  used  in  aeroelastic  analyses 
for  the  prediction  of  the  flutter  speed,  altitude,  and  frequency.  The 
variability  observed  in  the  test  data  is  then  reflected  through  a  similar 
variability  of  the  expected  flutter  characteristics.  It  also  implies  an 
uncertainty  on  where/when  this  phenomenon  might  occur  when 
using  components  that  have  not  been  tested. 

To  address  this  uncertainty  issue,  it  is  suggested  here  to  first 
proceed  with  a  structural  modeling  that  directly  accounts  for  its 
variability,  e.g.,  using  a  probabilistic  framework  as  will  be  done  here. 
Propagating  next  the  uncertainty  through  the  aeroelastic  com¬ 
putations,  typically  by  Monte  Carlo  methods,  will  then  provide  the 
desired  metrics  (e.g.,  standard  deviations,  98th  percentiles)  of  the 
flutter  boundary  variability.  In  this  regard,  note  that  the  consideration 
of  the  structural  uncertainty  can  be  carried  out  at  either  full  FEM  or 
modal  model  levels.  The  latter  is  adopted  here,  and  a  general 
framework  for  the  construction  of  a  stochastic  modal  model  is 
described  and  validated  using  experimental  results  on  the  Air  Force 
Institute  of  Technology  (AFIT)  joined  wing  [2]. 

II.  Stochastic  Modal  Model  Identification 

The  construction  of  a  stochastic  modal  model  of  a  structure 
requires  three  key  aspects:  1)  the  selection  of  a  basis,  including  its 
size  and  the  functions  that  compose  it;  2)  the  adoption  of  a 
representation  of  the  random  elements  (typically  the  elements  of  the 
modal  mass  and/or  stiffness  matrices)  of  the  modal  model;  and  3)  the 
estimation  of  the  deterministic  parameters  of  this  representation  from 
870  experimental  data. 
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Considering  first  the  basis  selection,  two  notable  alternatives  are 
possible.  The  first  mimics  the  corresponding  deterministic  analyses 
in  which  the  basis  functions  are  intimately  linked  to  the  structural 
model  (i.e.,  they  are  its  mode  shapes)  and  the  resulting  modal  mass 
and  stiffness  matrices  are  diagonal.  The  direct  extension  of  these 
concepts  to  the  case  of  a  stochastic  structure  leads  to  random  basis 
functions  with  a  modal  stiffness  (and  possibly  mass)  matrix  that 
remains  diagonal  with  random  elements.  This  representation  is 
certainly  acceptable  but  is  generally  very  challenging  to  establish 
because  of  the  need  to  characterize  the  joint  distribution  of  the 
elements  of  the  random  mode  shapes  guaranteeing  the  orthogonality 
and  normalization  conditions. 

A  simpler  approach  is  to  rely  on  a  fixed,  deterministic  basis  to 
which  will  then  be  associated  generally  full  stiffness  and/or  mass 
matrices.  This  strategy  is  often  employed  with  the  basis  functions 
being  the  mode  shapes  of  the  mean  model,  but  in  its  absence  here,  it  is 
proposed  to  adopt  the  mode  shapes  of  the  baseline  FEM,  as  they  span 
the  correct  spatial  domain  and  satisfy  the  proper  boundary 
conditions.  The  number  of  such  modes  to  include  in  the  basis  will  be 
estimated  from  the  experiments,  as  discussed  next. 

With  the  basis  fixed,  the  randomness  of  the  structure  is  localized  to 
the  elements  of  the  corresponding  stiffness  (and/or  mass)  matrix,  the 
statistical  distribution  of  which  must  then  be  characterized.  This  task 
is  performed  here  using  the  nonparametric  stochastic  modeling 
approach  (e.g.,  [3,4])  in  which  the  joint  probability  density  function 
of  the  random  (stiffness  and/or  mass)  matrix  elements  is  derived  to 
achieve  the  maximum  of  the  statistical  entropy  under  specified, 
physical  constraints.  The  nonparametric  approach  thus  replaces  the 
many  distribution-related  assumptions  usually  involved  in  the 
characterization  of  random  structural  models  by  a  single  hypothesis: 
the  maximization  of  the  entropy.  This  condition  implies  that  the 
distribution  of  the  uncertain  properties  will  have  a  heavy  tail,  i.e.,  will 
spread  as  widely  as  is  consistent  with  the  imposed  constraints.  Then, 
a  good  sampling  of  the  space  of  the  acceptable  values  of  these 
properties  will  be  achieved. 

The  constraints  imposed  on  the  maximization  of  the  entropy  stem 
from  physical  requirements  expected  of  all  stiffness  and  mass 
matrices:  symmetry,  positive  definiteness,  and  finiteness  of  the 
inverse  (assuming  that  there  are  no  rigid-body  modes).  It  is  further 
imposed  that  the  mean  matrix  be  equal  to  its  counteipart  for  the  mean 
model.  With  these  constraints,  it  was  shown  [3]  that  the  maximum  of 
the  entropy  is  maximized  when  the  random  matrices  B  to  be 
simulated  can  be  expressed  as 

B  =  LHHtLt  (1) 

where  L  is  any  decomposition  (e.g.,  a  Cholesky  decomposition)  of 
the  mean  matrix  B,  i.e.,  satisfying  B  =  LLl .  Furthermore,  H  denotes 
a  lower  triangular  random  matrix,  the  elements  of  which  are  all 
statistically  independent  of  each  other.  Moreover,  1)  the  diagonal 
elements  Hu  are  obtained  as  Hu  =  ^ Yti/ p, ,  where  Yu  is  Gamma 
distributed  with  parameter  (p(i)  —  l)/2  and  2)  the  offdiagonal 
elements  Hjh  i  Z,  are  normally  distributed  (Gaussian)  random 
variables  with  standard  deviation  a=l/s/2jl,  as  pictorially 
described  in  Fig.  1.  Note  further  that 

p(i)  =  n  —  i  +  2X  —  1  (2) 

and 


with  n  denoting  the  size  of  the  matrices  B. 

In  the  preceding  equations,  the  parameter  A,  >  0  is  the  free 
parameter  of  the  statistical  distribution  of  the  random  matrices  H  and 
B  and  can  be  evaluated  to  meet  any  given  information  about  their 
variability;  often,  the  standard  deviation  of  the  lowest  natural 
frequency  of  the  random  linear  system  will  be  specified  and  used  to 
obtain  X.  This  condition,  coupled  with  Eqs.  (1-3)  and  properties  1 
and  2,  provides  a  complete  scheme  for  the  generation  of  random 
symmetric  positive  definite  matrices  B. 


The  optimization-based  formulation  of  the  nonparametric 
approach  permits  the  consideration  of  different  or  additional 
constraints  as  required  by  the  problem  at  hand.  This  observation  has 
led  to  an  extension  [4]  of  the  preceding,  original  methodology  [3]  in 
which  the  standard  deviations  of  several  natural  frequencies  are 
imposed.  The  resulting  probabilistic  model  then  exhibits  more  than 
the  single  parameter  X  for  an  improved  matching  of  available  data. 
Furthermore,  a  computationally  efficient  simulation  strategy  of  the 
corresponding  random  matrices  B,  paralleling  but  different  from  the 
one  of  Eqs.  (1-3)  and  Fig.  1,  has  also  been  devised  [4], 

The  preceding  discussion  has  addressed  the  first  two  key  issues  in 
the  formulation  of  a  stochastic  modal  model  of  the  structures 
considered.  In  fact,  this  model  is  completely  characterized  by  the 
known  basis  functions  and  the  yet  unknown  mean  stiffness  matrix  K 
(assuming  that  the  uncertainty  is  in  stiffness,  and  thus  B  =  K)  and  the 
parameter  X  of  the  original  nonparametric  approach  (or  the 
parameters  A;  of  the  extended  formulation  [4]).  Note  that  the  neglect 
of  any  variability  in  the  mass  properties  leads  to  an  identity  modal 
mass  matrix  if  the  finite  element  (FE)  mode  shapes  have  been 
normalized  in  this  manner. 

The  last  issue  to  be  addressed  is  thus  the  estimation  of  the  matrix  K 
and  parameter  X  from  the  experimental  data.  In  this  context,  it  is 
assumed  here  that  M  natural  frequencies  (dm  and  mode  shapes  ipln 
of  a  series  of  R  similar  structures  have  been  measured:  m  = 
1 , 2, ....  M  and  r  =  1,2, ...  ,R.  Then,  before  proceeding  with  the 
estimation  of  K  and  X,  it  is  first  necessary  to  reduce  the  available  data 
to  the  modal  model  framework;  that  is,  the  mode  shapes  xpd  must  be 
projected  on  the  basis  of  N  (N  >  M)  modes  (pj  of  the  baseline  FEM. 
This  is  achieved  by  searching  for  the  best  representation  of  i/fd  in  the 
form 

(4) 
i=  i 

In  fact,  the  coefficients  ad  represent  the  mode  shape  m  of 
experiment  r  viewed  in  the  modal  space  spanned  by  the  FE 
modes  tf>  r 

If  all  degrees  of  freedom  (DOFs)  of  the  FEM  have  been  observed 
in  the  experiment,  the  deteimination  of  the  modal  coefficients  ad 
can  be  accomplished  by  premultiplying  Eq.  (4)  by  the  FE  mass 
matrix  and  using  the  orthogonality  properties  of  the  modes  (p  j .  If  not 
all  DOFs  have  been  observed,  as  is  expected  in  connection  with  most 
large  structures,  the  estimation  of  the  coefficients  ad  can  be  done  by 
a  least-squares  solution  of  the  linear  system  of  equations  [Eq.  (4)]. 
This  approach  will  be  referred  to  here  as  the  original  projection 
approach. 

The  overall  fit  of  the  approximation  of  Eq.  (4)  can  be  measured  by 
the  error 


zero  mean  Gaussian,  independent  of 
each  other  with  standard  dev.  a=l/^L 


and  10  (dev.  denotes  deviation). 
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where  the  subscript  N  emphasizes  the  dependence  of  the  error  on  N, 
the  number  of  FE  modes.  It  is  proposed  here  to  select  this  parameter 
to  be  the  lowest  value  for  which  En  is  below  a  specified  threshold 
value  related  to  the  noise  level  in  the  experiments. 

The  estimation  of  K  and  X  may  then  proceed  from  the  measure¬ 
ments  (dm  and  ad,  and  a  maximum  likelihood  approach  is  suggested 
here  [5],  as  demonstrated  in  [6-9]  forrelated  problems.  That  is,  K  and 
X  will  be  determined  to  yield  the  maximum  of  the  likelihood  function 

R 

L  =  Y\px(x(r)\K,K)  (6) 

r=  I 

where  px(x;K,X)  denotes  the  probability  density  function  of  the 
random  vector 


Z  —  (f2[  A[ |  . .  -  A ! } )  02  A2 1  . . .  A2(jv_ i)  ■  ■  •  ) 

Included  in  this  vector  are  the  M  random  natural  frequencies  Q,„, 

m=  1 . M,  of  measured  values  (dm  ;  and  the  random  mode  shapes 

coefficients  Amj,  m=  1 , ,M  and  j  =  1 . N  —  1,  of  measured 

values  dmj-  These  random  variables  are  derived  from  the  stochastic 
modal  model  with  random  stiffnesses  K  defined  by  Eqs.  (1)  and  (3) 
for  the  specified  K  and  X.  Note  that  only  the  first  N  —  1  component  of 
the  mode  shapes  are  included  in  X,  as  the  Nth  one  can  be  extracted 
from  the  normalization  condition 


N 

ATmAm  =  YJd1mj  =  1  (7) 

7=1 

given  that  the  modal  mass  matrix  is  identity.  Furthermore,  the  vector 
x{r)  denotes  the  realized  values  of  X  in  the  rth  experiment.  That  is, 

X  —  (&)  au  ■  •  ■  tt>2  “21  ■  •  ■  a2(JV-l)  ■  ■  •  ) 

The  estimation  of  px(x;  X,  X)  would  proceed  as  follows.  For  a 
given  mean  stiffness  matrix  K  and  uncertainty  parameter  X,  a  large 
number  of  realizations  of  the  random  modal  stiffness  matrix  K  would 
be  generated  according  to  Eqs.  (1—3).  Then,  for  each  one  of  them,  the 
natural  frequencies  and  mode  shapes  of  the  corresponding  dynamical 
system  would  be  obtained.  Then,  a  statistical  analysis  of  this 
population  of  frequencies  and  mode  shapes  would  lead  to  the 
estimation  of  px(x;  K,  A).  Given  the  expected  high  dimensionality  of 
X,  with  p  =  N  x  M  components,  the  determination  of  this  exact 
probability  density  function  may  represent  a  daunting  task  requiring 
an  excessively  large  number  of  simulations  to  be  carried  out.  An 
alternate  approach  is  to  approximate  px(x\  K ,  X)  by  its  Gaussian 
counterpart  (asymptotically  correct  as  the  uncertainty  level  goes  to 
zero;  see  [10]): 


Px(x\  K,  X) 


1 

(2ny'2^fei(Kdd) 


III.  Test  Article:  Air  Force  Institute  of  Technology 
Two  Meter  Joined  Wing 

The  development  of  persistent  intelligence,  surveillance,  and 
reconnaissance  (ISR)  platforms  is  a  key  area  of  emphasis  for  the 
foreseeable  future.  Current  airborne  ISR  assets  can  only  provide 
limited  persistent  coverage,  and  their  current  sensors  can  only  collect 
data  over  a  limited  area  when  compared  with  radar.  The  Boeing 
concept  includes  a  joined-wing  design,  as  shown  in  Fig.  2  [14].  The 
joined-wing  design  concept  allows  an  aircraft  to  have  a  large 
wingspan  due  to  the  inclusion  of  an  aftwing,  which  provides  the 
additional  required  stiffness.  The  long  wingspan  provides  the  ability 
to  house  a  large  radar  in  order  to  generate  long  wavelengths 
necessary  to  penetrate  foliage.  These  radar  arrays  will  have  to  be  fully 
integrated  into  the  remotely  piloted  vehicle  wing  structure  in  order  to 
keep  size  and  weight  minimal.  This  concept  of  combining  the  radar 
array  and  structure  is  referred  to  as  the  conformal  load-bearing 
antenna  structure. 

To  achieve  high  lift  and  low  drag,  the  wings  will  have  a  very  high 
aspect  ratio.  This  concept  vehicle’s  long  endurance  requirement 
dictates  a  high  fuel  load,  a  low  structural  weight  fraction,  and  a  high 
aspect  ratio,  and  it  will  likely  result  in  wings  that  are  flexible, 
resulting  in  large  deflections  under  load.  The  aftwing  is  present  to 
improve  structural  efficiency  in  achieving  larger  wingspans  without 
significantly  increasing  the  structural  weight  required  to  withstand 
high-lift  loads.  It  is  expected  that  this  large  wingspan  will  have  a  high 
geometric  nonlinearity  with  respect  to  external  loads  [15,16]. 

Wolkovich  presented  the  joined-wing  concept  in  1986  [17].  In  his 
paper,  Wolkovich  demonstrates  two  important  advantages  of  the 
joined-wing  concept:  light  weight  and  high  stiffness.  Early  studies 
showed  that  bending  due  to  a  vertical  tip  load  occurs  in  the  direction 
perpendicular  to  the  plane  created  by  the  axes  of  the  fore-  and 
aftwings  [18].  Figure  3  depicts  the  plane  created  by  the  axes  of  the 
fore-  and  aftwings.  For  long,  flexible  joined  wings  under  a  large  load, 
the  applied  tip  load  will  cause  bend-twist  coupling  in  the  aftwing. 
However,  Patil  noted  that  the  joined-wing  design  is  much  stiffer  than 
single-wing  configurations,  and  in  his  tests,  he  also  found  that  the 
joined  wing  can  be  made  stiff  enough  to  only  show  negligible 
nonlinear  effects  [19]. 

Since  the  joined  wing’s  inception,  there  has  been  a  wide  variety  of 
research  efforts  on  the  joined-wing  concept.  Snyder  et  al.  analyzed  a 
twin-boom  fuselage  sensor-craft  configuration  that  included 
aeroelastic  analysis  using  FE  modeling  and  computational  fluid 
dynamics  [20],  Blair  et  al.  investigated  aspects  of  the  joined-wing 
structure  while  accounting  for  critical  loads  throughout  the  mission 
envelope  [21].  Smallwood  investigated  the  joined  wing  and 
predicted  a  worst-case  gust  load  response  resulting  in  about  an  8- 
9  deg  pointing  error  [22] .  Roberts  computed  predictions  to  a  gust  load 
response  while  accounting  for  geometric  nonlinearities  [23].  Bond 
and  Canfield  successfully  created  an  all-aluminum  2  m  semispan  test 
article  on  which  they  applied  static  follower  forces  and  collected 
nonlinear  deformations  [24].  They  were  also  successful  in  reaching 
agreement  between  their  measured  displacements  and  strains  and 
Nastran  nonlinear  analytical  results.  Based  on  lessons  learned  from 
Bond  and  Canfield’s  work,  a  2  m  joined-wing  test  article  was 
designed  and  built  at  AFIT,  which  is  discussed  next. 

The  primary  test  article  discussed  in  this  research  effort  is  a  2  m 
joined  wing,  as  shown  in  Fig.  4,  which  has  a  winglike  cross  section. 
The  fore-  and  aftwings  are  solid  aluminum  with  cross-section 


x  exp 


1 

2 


(x  -  fix)TKx'x(x  -  px) 


(9) 


where  fix(K ,  X)  and  KXX(K,  X)Kxx(K,  X)  denote  the  mean  vector 
and  covariance  matrix  of  the  random  vector  X  determined  from  the 
simulations  described  previously  for  the  given  mean  stiffness  matrix 
K  and  uncertainty  parameter  X.  The  Gaussian  assumption  of  Eq.  (9) 
leads  to  a  computationally  very  efficient  approach  but  may  not 
always  provide  the  desired  close  fit  of  the  probability  density 
function.  In  such  situations,  the  more  refined  kernel  density 
estimation  method  [11-13]  could  be  used  instead  of  Eq.  (9)  to  build 
the  likelihood  function. 


Fig.  2  Boeing  sensor-craft  concept. 
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Fig.  4  2  m  joined-wing  test  article. 


dimensions  of  1  by  8  in.  for  the  forewing  and  6  in.  for  the  aft  wing.  All 
attachment  fixtures  are  made  of  hardened  4130  tool  steel  to  maximize 
safety  and  minimize  motion  at  the  wing  root  attachment  points 
attached  to  the  concrete  floor.  Besides  approximating  the  geometric 
layout  of  the  joined-wing  configuration,  the  original  primary  design 
focus  for  static  testing  was  the  inclusion  of  a  nonlinear  bend-twist 
couple  without  permanent  deformation  of  the  test  article. 

For  modal  testing,  a  Polytec  PS  V-400-3D  scanning  laser  Doppler 
vibrometer  (LDV)  was  used  to  collect  operational  deflection  shapes, 
natural  frequencies,  and  damping  ratios. 

IV.  Finite  Element  Model 

Three  distinct  FEMs,  each  composed  of  nearly  all  beam,  plate,  and 
solid  elements,  of  the  2  m  joined-wing  test  article  were  originally 
created.  The  predominantly  plate  element  model  has  almost 
60,000  DOFs  (see  Fig.  5)  and  yielded  the  best  results  in  terms  of 
minimal  model  complexity  and  agreement  between  the  measured 
and  predicted  static  displacements.  The  hardened  steel  mounting  and 
joint  structures  are  modeled  as  well,  while  the  concrete  is  assumed  to 
be  a  rigid  constraint.  Connections  between  all  plates  are  represented 
with  rigid-link  and  6-DOF  spring  elements.  The  6-DOF  spring 
elements  allow  a  less  than  rigid  connection  between  components. 


The  values  of  these  stiffnesses  were  initially  selected  [2]  to  match 
experimentally  measured  static  deflections.  The  stiffnesses  of  the 
spring  elements  of  the  connection  between  the  wings  and  the  bottom 
mounting  structures  were  tuned  a  second  time1  to  provide  a  closer  fit 
of  the  first  natural  frequency  with  its  measured  counterpart.  This 
model,  without  any  further  tuning,  was  considered  as  the  baseline 
model  discussed  in  the  previous  section  and  led  to  natural 
frequencies  shown  in  Table  1  and  mode  shapes  as  depicted  in  Fig.  6 
[projected  on  the  laser  direction;  see  discussion  of  Eq.  (11)]. 

V.  Experimental  Procedure 

Figure  5  shows  the  252  measurement  locations  as  dots  on  the 
structure.  These  points  were  scanned  using  the  Polytec  PSV-400-3D 
scanning  LDV.  The  joined-wing  test  article  was  excited  with  an 
autoping  hammer  with  a  force  sensor  mounted  to  the  tip  of  the 
stinger.  The  excitation  location  remained  fixed  for  all  tests. 
Frequency  response  functions  (FRFs)  were  computed  for  all 
measurement  locations  after  10  averages.  The  252  FRFs  collected 
were  spectrally  sieved  around  the  resonance  peaks,  and  a  modal 
parameter  estimation  program  (MEscope)  was  used  to  extract 


^Personal  communication  with  B.  Grier  and  R.  Figliola,  2010. 
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Fig.  5  2  m  joined-wing  FEM.  Dots  represent  the  252  measurement 
locations  that  are  coincident  with  the  FEM  nodes. 


accurate  estimates  of  the  mode  shapes,  natural  frequencies,  and 
modal  damping  factors.  Note  that  the  extracted  mode  shapes  were 
complex,  although  very  close  to  purely  imaginary  for  the  points  with 
significant  motion.  Thus,  the  imaginary  parts  and  absolute  values  of 
the  modes  could  both  be  used  as  estimates  of  the  mode  shapes.  A 
third  set  of  such  estimates  was  obtained  using  the  Niedbal  approach 
[25],  in  which  the  complex-valued  modal  matrix  \PC  (the  columns  of 
which  are  the  estimated  complex  modes)  is  transformed  into  a  real¬ 
valued  one  'f'/v  according  to 

*jv  =  ReWc)  +  rm(yc)[Re(yc)TRe(yc)]-'Re(yc)TIm(yc) 

(10) 


A.  Structural  Population 

A  key  in  the  present  effort  is  the  availability  of  test  data  on  different 
structures  for  the  estimation  of  the  level  of  structural  uncertainty.  As  a 
first  validation  of  the  proposed  concepts  and  methods,  it  was  decided 
to  simulate  the  variability  in  structures  by  disassembling  and 
reassembling  the  wings.  In  this  process,  the  22  bolts  of  the  top  plate  of 
the  top  mounting  structure  were  first  removed,  releasing  the  two 
wings,  which  were  separated  from  each  other  by  6  mm.  Then,  these 
two  wings  were  brought  back  together,  the  top  plate  was  reinserted, 
and  the  22  bolts  were  refastened  to  the  specified  torque. 

The  testing  described  in  the  previous  subsection  was  first 
performed  on  the  wing  in  the  configuration  that  it  was  left  in  after  the 
testing  reported  in  [2],  To  obtain  an  assessment  of  the  variability  in 
the  measurement  and  data  analysis  process,  three  repeat  tests  were 
carried  out  on  that  structure  (referred  to  as  structure  1  in  the 
remainder  of  this  paper)  without  any  disassembly.  The  extracted 
natural  frequencies  are  given  in  Table  1  (columns  labeled  Repeat  1 , 
Repeat  2,  and  Repeat  3).  The  disassembly  and  reassembly  processes 
were  then  carried  out  three  times  to  obtain  structures  2,  3,  and  4, 


which  were  each  tested  as  described  previously  but  only  once;  see 
Table  1  for  the  corresponding  natural  frequencies  and  Fig.  6  for  the 
mode  shapes  of  structure  4. 

B.  Modal  Projection  on  the  Finite  Element  Basis 

The  first  step  in  the  estimation  of  the  mean  matrix  K  and 
uncertainty  parameter  X  is  the  projection  of  the  experimental  mode 
shapes  on  the  FE  basis  [see  Eq.  (4)]  and  the  selection  of  the  number  N 
of  such  modes  to  be  included  in  the  basis.  In  this  regard,  it  must  first 
be  noted  that  a  single  laser  vibrometer  was  used  in  the  experimental 
setup;  thus,  the  mode  shapes  identified  correspond,  in  fact,  to  the 
three-dimensional  modal  deflections  projected  along  the  direction  of 
the  laser  beam,  from  the  collimator  to  the  measurement  point.  To 
maintain  consistency  in  the  estimation  of  the  coefficients  the 
basis  <j)j  was  then  constructed  using  a  similar  projection  of  the  three- 
dimensional  FE  modes  on  the  direction  of  the  local  direction  of  the 
laser  beam.  That  is, 


(Sj-SL) 

Ik- II 


(11) 


where  (f>jj  and  fM  are  the  components  of  the  modal  basis  function 
and  FE  mode  shape  j  corresponding  to  measurement  point  i  of 
coordinate  s,  .  Furthermore,  sL  denotes  the  position  vector  of  the  laser 
collimator.  The  computations  of  Eq.  (11)  were  carried  out  for  FE 
modes  1, 2, 3, 4, 5, 7,  etc.  Mode  6  was  not  included  in  this  list  because 
it  is  predominantly  in-plane  (at  the  contrary  of  the  others 
predominantly  transverse)  and  because  it  was  not  observed 
experimentally  in  the  band  of  interest  (i.e.,  below  100  Hz). 

The  projection  of  the  experimental  mode  shapes  on  the  basis  (j>  ■ 
(see  Fig.  6  for  structure  4)  was  next  carried  out  for  repeat  3  of 
structure  1  and  the  measurements  of  structures  2-4  for  an  increasing 
number  N  of  transverse-dominant  basis  functions,  and  the  overall 
error  [see  Eq.  (5)]  was  determined  in  each  case.  The  evolution  of  this 
error  with  N  is  shown  in  Fig.  7  and  suggests  that  the  use  of  N  =  6  FE 
modes  (modes  1,  2,  3, 4,  5,  and  7)  allows  the  error  to  reach  its  noise- 
related  floor.  Note  that,  in  Fig.  7,  experimental  mode  5  was  not  used 
in  the  computation  of  the  error  because  it  is  significantly  more  noisy 
than  the  other  modes  owing  to  the  accidental  coincidence  of  its 
natural  frequency  with  the  60  Hz  electric  power  frequency.  The 
inclusion  of  this  mode  in  the  error  raises  significantly  the  noise  floor 
level  and  renders  more  difficult  the  capturing  of  the  convergence  of 
the  error.  Note  finally  that  the  results  presented  in  Fig.  7  were 
obtained  with  the  experimental  mode  shapes  estimated  by  all  three 
methods:  imaginary  part,  absolute  of  the  complex  mode,  and  Niedbal 
transformation  [25] ;  see  Eq.  ( 1 0).  The  first  two  provide  errors  that  are 
almost  exactly  equal,  while  the  latter  gives  a  very  slightly  larger 
value.  In  all  ensuing  discussions,  the  imaginary  part  of  the  complex 
mode  was  used. 

With  the  number  of  basis  functions  and  the  coefficients  a„j 
estimated,  it  is  possible  to  form  reconstructed  mode  shapes,  i.e.,  the 
projections  of  the  experimental  ones  on  the  basis  as  expressed  by  the 
right-hand  side  of  Eq.  (4).  These  modes  are  also  shown  in  Fig.  6  for 
the  experimental  data  of  structure  4.  The  modal  assurance  criterion 
(MAC)  values  of  the  experimental  data  (imaginary  part  of  the 
complex  modes)  with  the  corresponding  basis  function  (FE  mode 
projected  on  the  laser  direction)  and  with  the  reconstructed  mode 


Table  1  Natural  frequencies  (Hz)  obtained  from  the  FEM  and  experimentally  measured" 


FEM 

Repeat  1 

Repeat  2 

Repeat  3 

Struct.  2 

Structure  3 

Structure  4 

1 

8.30 

8.33 

8.31 

8.30 

8.26 

8.29 

8.27 

2 

16.09 

15.33 

15.38 

15.37 

14.95 

15.03 

15.00 

3 

19.50 

20.61 

20.58 

20.58 

20.54 

20.58 

20.58 

4 

33.64 

38.08 

38.07 

38.05 

37.70 

37.74 

37.07 

5 

55.45 

60.04 

59.97 

59.96 

59.98 

59.92 

59.88 

6 

73.93 

>100 

>100 

>100 

>100 

>100 

>100 

7 

75.51 

88.37 

88.55 

88.50 

86.03 

86.22 

80.50 

“Repeats  1,  2,  and  3  denote  three  tests  of  structure  1. 
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Fig.  6  Experimental  mode  shapes  (imaginary  part  estimate),  laser  projected  FE  mode  shapes,  and  experimental  mode  shapes  projected  in  the  basis. 


shapes  are  also  presented  on  Fig.  6.  The  high  values  (greater  than  or 
equal  to  0.972)  of  the  former  MAC  (between  experiment  and  basis 
functions)  demonstrate  the  good  accuracy  of  the  FEM  (even  though 
there  are  notable  discrepancies  on  natural  frequencies;  see  Table  1) 
and  the  low  level  of  noise  in  the  laser  data.  After  reconstruction,  the 
MAC  values  are  even  higher  (greater  than  or  equal  to  0.993),  clearly 
confirming  the  adequacy  of  the  basis  for  the  representation  of  the 
experimental  mode  shapes  and  the  presence  of  a  slight  coupling 
between  the  basis  functions. 


This  coupling  can  also  be  observed  in  the  projection  coefficients 
alnj,  which  are  very  close  to  1  for  m  =  j  and  rather  small  for  m  ^  j; 
see  Table  2  for  structure  1,  repeat  3  data,  which  is  rather  typical.  A 
pictorial  representation  of  these  coefficients  is  shown  in  Fig.  8.  Note 
from  this  figure  and  from  Table  2  that  the  column  associated  with 
experimental  mode  5  (of  frequency  very  close  to  60  Hz)  involves 


larger  coefficients  cYY  for  m  ^  j ,  owing  probably  to  the  larger  noise 
level  on  this  mode. 


Approved  for  public  release;  distribution  unlimited. 


6 


876 


AVALOS  ETAL. 


Fig.  7  Original  projection  error  of  Eq.  (5)  vs  number  of  transverse- 
dominant  basis  functions.  Experimental  modes  5  and  6  not  used.  Mode 
shape  selected  from  imaginary  part,  absolute  value  of  the  complex  mode, 
or  through  the  Niedbal  transformation. 


VI.  Mean  Model  and  Uncertainty 
Parameter  Estimation 

The  estimation  of  the  mean  stiffness  matrix  K  and  uncertainty 
parameter  X  from  the  measured  natural  frequencies  and  projection 
coefficients  according  to  the  maximum  likelihood  approach  is 
estimated  with  a  base  of  six  FEM  modes. 

A.  Optimization  Algorithm  Specificities 

As  can  be  seen  from  Eqs.  (6-9),  the  evaluation  of  the  likelihood 
function  requires  the  determination  of  the  mean  vector  fLx(K,k)  and 
the  covariance  matrix  KXX(K,  X).  Yet,  neither  one  can  be  determined 
in  closed-form  expression;  thus,  it  is  necessary  to  proceed  with  a 
Monte  Carlo  simulation.  Let  S  be  the  number  of  random  stiffness 
matrices  simulated  for  a  particular  mean  matrix  K  and  uncertainty 
parameter  X.  For  each  of  these  matrices,  an  eigenvalue  computation  is 
carried  out  and  yields  a  vector  X  =  Y(I)  composed  of  the  natural 
frequencies  and  mode  shape  values.  Then,  the  mean  fix  of  this  vector 
can  be  estimated  from  the  population  of  size  S  as 

(12) 

5=1 

The  covariance  matrix  is  similarly  estimated  from  the  samples 
Z(I)  as 

Kxx  (X<S>  -  ^)T  (13> 

5=1 

Because  of  the  finiteness  of  the  number  of  samples  5,  Eqs.  (12)  and 
(13)  only  provide  estimates  of  the  true  mean  and  covariance  matrix. 
Furthermore,  these  estimates  are  random;  they  vary  from  one 
population  of  5  samples  to  another  exhibiting  a  standard  deviation 
that  is  typically  of  the  order  of  1  /  s/~S. 


Fig.  8  Projection  coefficients  a^j  for  structure  1,  repeat  3. 

Thus,  if  a  perturbation  A  in  either  K  or  X  is  performed, 
corresponding  changes  in  fix(K,  X)  and  KXX(K,  X)  are  obtained  that 
are  composed  of  two  parts.  The  first  component  is  the  change  in  the 
true  value  of  the  mean  vector  and  covariance  matrix  that  is  implied  by 
A.  The  second  part,  however,  is  a  random  term  that  results  from  the 
finiteness  of  S  and  the  approximate  nature  of  Eqs.  (12)  and  (13).  This 
latter  term  renders  difficult  the  accurate  evaluation  of  gradients, 
especially  when  they  are  close  to  zero,  unless  the  sample  size  S  is 
taken  to  be  prohibitively  large. 

These  observations  demonstrate  that  the  numerical  optimization 
of  the  likelihood  function  should  be  accomplished  by  relying  on 
methods  that  do  not  need  gradient  evaluation,  as  is  necessary  when 
the  objective  function  is  not  smooth.  In  this  effort,  pattern  search  and 
simulated  annealing  methods  were  employed  for  finding  the 
maximum  of  the  likelihood  function.  In  fact,  simulated  annealing 
was  found  to  perform  better  in  finding  the  best  (global)  optimum  than 
pattern  search  did. 

A  second  specificity  of  the  current  optimization  problem  concerns 
the  mean  stiffness  matrix  K.  Indeed,  it  is  necessary  that  this  matrix 
remains  positive  definite  and  symmetric.  While  the  latter  property  is 
easily  built  in,  the  former  is  not.  Imposing  the  positiveness  of  the 
eigenvalues  of  K  leads  to  algorithmic  difficulties,  and  an  alternative 
formulation  is  desired.  The  constrained  optimization  can  be  recast 
into  an  unconstrained  one  by  varying,  not  the  components  of  K  but 
rather  the  components  of  the  lower  triangular  matrix  LK  corre¬ 
sponding  to  the  Cholesky  decomposition  of  K,  i.e.,  with 

K  =  LkLtk  (14) 

In  fact,  this  decomposition  is  necessary  in  the  nonparametric 
modeling  effort.  Note  that  the  results  presented  herein  were  obtained 
with  a  slightly  different  formulation  in  which  the  mean  stiffness 
matrix  is  expressed  as  would  be  obtained  from  its  eigenvalues  and 
eigenvectors,  i.e.,  as 

K  =  Q-  A  - Qt  (15) 

where  A  is  diagonal  with  positive  elements. 


Table  2  Projection  coefficients  for  structure  1,  repeat  3 


Experimental  modes 

Basis  function 

i 

2 

3  4 

5 

7 

1 

0.9989 

-0.1340 

0.1159 

0.0310 

0.2200 

0.0455 

2 

-0.0386 

-0.9861 

0.0289 

-0.0889 

-0.0750 

-0.1695 

3 

0.0236 

0.0358 

-0.9911 

0.0758 

0.1475 

0.0576 

4 

-0.0093 

0.0793 

-0.0245 

-0.9899 

0.1202 

0.2432 

5 

0.0065 

0.0266 

0.0467 

0.0343 

-0.9421 

0.0001 

6 

0.0023 

0.0304 

0.0016 

-0.0035 

0.0173 

-0.9493 

7 

-0.0047 

0.0074 

-0.0266 

0.0011 

-0.1180 

-0.0227 

8 

0.0050 

0.0075 

0.0020 

-0.0113 

0.0547 

0.0234 

9 

0.0015 

0.0117 

-0.0022 

0.0636 

0.0707 

0.0574 

10 

-0.0015 

0.0135 

-0.0021 

-0.0116 

0.0061 

-0.0362 
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B.  Normalization 

Numerical  difficulties/instabilities  were  encountered  in  the 
optimization  process  when  the  uncertainty  level  was  small  (i.e.,  large 
X)  and  which  were  traced  to  the  small  values  of  the  elements  of  the 
covariance  matrix.  This  issue  can  be  eliminated  by  appropriately 
renormalizing  the  variables.  Specifically,  proceed  first  with  the 
Cholesky  decomposition  of  the  covariance  matrix  as 

Kxx  =  (16) 

Next,  introduce  the  new  random  variables: 

Y  =  V~l(X  —  fix)  (17) 


Then,  proceeding  with  the  change  of  random  variables  leads  to 


p$(y\K,X) 


1 

(2  n)p/2 


exp 


(18) 


Fig.  9  Maximum  value  of  the  log-likelihood  as  a  function  of  uncertainty 
level  (3),  original  method. 


since  the  covariance  matrix  KYy  is  readily  shown  to  be  identity. 

C.  Likelihood  Function  Versus  Log-Likelihood  Function 

Rewriting  the  likelihood  function  in  terms  of  the  variables  Y  leads 
to 


L  =  f\Pr(y^K,X)^=f\Pr(y^K,X) 


■s/det(Kxx) 
(19) 


From  the  orthogonality  property,  one  concludes  that  [AM]  1  = 
Awr  or 

[I  +  AA<r)]-'  =  [I  +  A  Awf  =  /  +  AA<r,r  (24) 

Using  a  Neumann  expansion  for  small  A Air\  i.e.,  such  that  all  of 
its  eigenvalues  are  significantly  less  than  1  in  magnitude,  one  obtains 

[7+  AAw]-‘  «/-AA«  (25) 

Comparing,  finally,  Eqs.  (24)  and  (25)  leads  to 


Occasionally,  small  values  of  the  probability  density  function  are 
obtained  that  can  result  in  numerical  difficulties  (underflow).  This 
issue  is  resolved  by  proceeding  with  the  maximization  of  the  log- 
likelihood  function  log(L),  which  is  obtained  as 

^)-H  001 

D.  Orthogonality  of  Modes 

The  mode  shapes  simulated  on  the  basis  of  the  nonparametric 
model  correspond  to  a  symmetric  eigenvalue  problem  and  are  thus 
orthogonal  with  respect  to  both  mass  and  stiffness  matrices.  In  the 
present  effort,  the  modal  basis  was  formed  from  mass  normalized 
mode  shapes;  thus,  the  mean  (constant)  mass  matrix  of  the  reduced- 
order  model  is  identity.  Thus,  all  simulated  values  a satisfy  the 
orthogonality  condition 


A  A^T  =  -A  A«  (26) 

which  demonstrates  that  AA<r)  is  skew  symmetric  because  of  the 
orthogonality  condition  when  the  uncertainty  level  is  small.  While 
this  property  does  not  hold  for  larger  uncertainty  levels,  the  preceding 
derivation  nevertheless  demonstrates  that  there  is  a  strong  relation 
between  the  upper  and  lower  triangular  elements  of  A A(r) ;  thus,  only 
one  of  these  two  groups  should  be  considered  independent  variables. 

In  this  light,  the  likelihood  function  was  not  built  on  the 
distribution  of  the  entire  set  of  coefficients  ctj8],  j  =  1 , ,N  —  1, 

skew-symmetric 

-0.0396  0 

-0.0073  0.0132  0 

^mJ  ~  -0.0379  0.0633  -0.0636  0 

(Jf08f2>Ca0855^  0.0621  -0.0496  0 

'jjf0&7D>  0.0466  0.0034  -0.0697  -0.0225  0 


A«rA(r)  _  / 


(21) 


a) 


where  A(r)  is  the  matrix  of  the  coefficients  qA].  Such  a  property  is  in 
general  not  satisfied  by  the  experimental  measurements  because  of 
noise  and  differences  between  the  real  structure  and  the  compu¬ 
tational  model.  This  lack  of  consistency  between  the  properties  of  the 
simulated  and  experimental  modes  was  found  to  be  a  significant 
problem  in  the  reliable  estimation  of  the  uncertainty  parameter  X  or 
equivalently  of  S  defined  as 


0 

0.0689 


P  mi  = 


0 


-0.0288  0.0005  0 

-0.0552  0.0296  -0.0585 
0.0270  0.0848 

3jj)59T>  0.0235 


skew-symmetric 


0 

-0.014  0 

-0.0554  0.0028  0 


b) 


s2  = 


n  +  1 
n  +  2X  —  1 


(22) 


The  orthogonality  condition  existing  on  the  simulated  mode 
shapes  implies,  in  fact,  that  the  coefficients  cannot  be  indepen¬ 
dent  of  each  other.  In  fact,  note  that  the  matrices  Aw  are  very  close  to 
the  identity  matrix  for  a  small  level  of  uncertainty;  that  is,  they  can  be 
expressed  as 


0 


0.084 

0.0429 

0 

0.0345 

0 

symmetric 

0.0245 

0.0131 

0.0076 

0 

0.0529 

0.0729 

0.0328 

0.0692 

0 

0.0439 

0.0432 

0.0189 

0.0541 

0.017  0 

Aw  =  /  +  AA(r) 


Fig.  10  Representations  of  a-b)  values  of  pmj  obtained  in  the  first  and 
(23)  third  tests,  and  c)  standard  deviations  of  the  random  coefficients  . 
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m  =  1, ...  ,N,  as  suggested  by  Eqs.  (6-8)  but  rather  of  the 
combinations 

Pmj  =  -  “>!l/2  for  j  =  1 . m  -  1 ,  m  =  2, . . . ,  N  (27) 

Using  this  set  of  variables  and  the  orthogonality  conditions  of 
Eq.  (21),  the  entire  matrix  A(r)  can  be  uniquely  built.  Thus,  the  use  of 
the  set  [qA*  —  qA]/2  takes  automatically  into  account  the 
orthogonality  constraints,  which  thus  need  not  be  considered  any 
further. 

E.  Effect  of  Noise  and  Modal  Coupling 


The  preceding  optimization  process  was  applied  to  the  joined- 
wing  data.  The  initial  condition  used  for  the  mean  stiffness  matrix  K 
was  selected  as  the  matrix  for  which  the  eigenvalues  are  the  squared 
frequencies  average  over  the  tests  considered  and  with  mode 
shapes  similarly  averaged  from  the  matrices  A.  Furthermore,  the 
optimization  process  was  carried  out  for  a  set  of  fixed  values  of  the 
uncertainty  parameter  (iV  or  6)  and  the  resulting  maximum  value  of 
the  log-likelihood  function  was  determined  for  each  uncertainty 
level.  Collecting  these  results  led  to  the  plot  of  Fig.  9,  which  shows 
that  the  highest  value  of  the  likelihood  function  is  obtained  for  a 
rather  large  uncertainty  level,  i.e.,  S  =  0.6. 

This  optimum  value  of  the  uncertainty  level  appears  significantly 
larger  than  would  have  initially  been  expected  from  the  natural 
frequencies.  The  spread  in  these  values  seems  consistent  with  a  value 
of  5  =  0.1  or  less.  This  contradiction  suggests  that  the  high 
uncertainty  level  estimated  results  from  some  peculiarity  of  the  test 
data  mode  shapes.  To  obtain  a  better  perspective  on  this  situation, 
shown  in  Fig.  10  are  the  values  of  j8mJ  obtained  in  the  first  and  third 
tests,  as  representative  examples,  as  well  as  the  standard  deviation  of 
these  values  obtained  from  the  nonparametric  approach  with  an 
uncertainty  level  of  S  =  0.3. 

An  inspection  of  the  test  values  of  Figs.  1 0a  and  1  Ob  in  comparison 
with  the  standard  deviations  of  Fig.  10c  demonstrate  that  there  are  a 
few  outliers  in  the  data,  i.e.,  values  that  exceed  the  three  standard 
deviations  threshold.  Clearly,  the  maximum  likelihood  method  will 
be  biased  by  these  values  and  will  tend  to  increase  the  level  of 
uncertainty  in  the  system  until  these  outliers  have  a  higher  probability 
level.  An  equivalent  perspective  can  be  drawn  from  the  plots  of 
Fig.  11,  which  provide  a  series  of  two-dimensional  joint  probability 
density  functions  of  the  natural  frequencies  and  values. 

Plotted  on  these  figures  are  1)  the  test  points  (as  large  stars), 
2)  one-,  two-,  and  three-standard-deviation  contour  lines  [which  are 
ellipses  according  to  the  approximate  Gaussian  distribution  of 
Eq.  (8);  shown  as  lines],  and  3)  a  series  of  simulated  values  of  the 
frequencies  and  /ijd  values  (as  small  circles).  The  presence  of  the 
outliers  is  clearly  shown  as  well  on  these  plots. 

Among  the  outliers,  one  recognizes  particularly  the  terms  (5,1), 
(5,2),  (6,1),  and  (6,2),  and  these  components  surface  as  outliers  in 
particular  because  of  the  small  values  of  the  corresponding  standard 
deviations.  In  fact,  an  inspection  of  Fig.  10c  reveals  that  the  standard 
deviations  are  at  a  maximum  between  modes  of  consecutive  natural 
frequencies  and  steadily  decrease  as  the  separation  between  natural 
frequencies  increases. 

Why  then  are  some  of  these  values  so  large?  Those  associated 
with  mode  5  can  easily  be  dismissed  as  noise  related.  Indeed,  the 
natural  frequency  of  this  mode  is  almost  exactly  60  Hz,  and  a  strong 
presence  of  the  power  frequency  in  the  response  was  observed  at  the 
time  of  the  test.  Accordingly,  a  good  practice  would  be  not  to  involve 
this  mode  in  the  computations.  This  argument  then  leaves  the  (6,1) 
and  (6,2)  components.  The  former  component  may  exhibit  a  trend, 
being  large  in  the  last  three  tests  and  of  the  same  sign.  The  latter  one, 
however,  appears  only  in  the  last  test.  Since  the  values  of  these 
coefficients  /3mj  are  not  exceedingly  large,  it  is  indeed  possible  that 
noise  is  still  responsible  for  the  issue. 

To  address  the  presence  of  the  outliers,  two  different  efforts  were 
carried  out  in  this  investigation.  First,  the  modeling  procedure  was 


Fig.  12  Maximum  value  of  the  log-likelihood  as  function  of  the 
uncertainty  level  (6):  original  method  without  mode  5. 


repeated  without  mode  5  to  eliminate  the  effect  of  the  noise  implied 
by  the  coincidence  of  its  natural  frequency  with  the  power  frequency. 
This  process  did  not  lead  to  a  significant  change  in  the  optimum  value 
of  5,  most  likely  because  of  the  remaining  presence  of  the  (6,1)  and 
(6,2)  components;  see  Fig.  12. 

A  second  effort  focused  on  reformulating  the  projection  of  the  test 
mode  shapes  on  their  FE  counterparts  to  minimize  the  effect  of  noise. 
This  revised  formulation  is  described  next. 


F.  Revised  Projection  Formulation 

The  determination  of  the  projection  coefficients  amj  has  been 
performed  up  to  this  point  to  minimize  the  representation  error 


F(m)  _ 

^rep  — 


(28) 


for  each  mode  m  of  test  r.  If  noise  is  present  in  the  mode  shape 
measurements,  it  will  be  modeled  as  part  of  the  test  mode  and  will  be 
present  in  the  ensuing  computations,  potentially  biasing  the  maxi¬ 
mum  likelihood  identification  strategy  toward  a  higher  uncertainty 
level. 

In  this  light,  it  would  be  desirable  to  filter  out  the  noise,  i.e.,  to  fine- 
tune  the  estimation  of  the  coefficients  toward  their  use  in  the 
nonparametric  modeling  effort.  As  in  the  original  projection-based 
determination  of  the  coefficients  amj,  this  process  will  be 
accomplished  mode  per  mode  and  test  per  test.  For  a  given  mean 
model  K  and  uncertainty  level  X ,  the  fitness  of  a  set  of  coefficients  amj 
can  be  assessed  through  the  discrepancy  measure 


&  =  E 


(“mi  -  “my)2 


(29) 


where  amj  and  amj  are  the  mean  and  standard  deviations  of  the 
random  mode  shapes  AmJ  obtained  with  K  and  X.  The  larger  F^]  is, 
the  less  fit  the  coefficients  amj  are  for  nonparametric  modeling  or, 
equivalently,  the  larger  the  corresponding  uncertainty  level  will  be. 
Accordingly,  it  is  desired  to  minimize  Fm  ■  However,  this 


Fig.  13  Iterative  process  to  obtain  the  projection  coefficients,  revised 
formulation  (MC  denotes  Monte  Carlo,  and  NP  denotes  nonparametric 
model). 
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Optimum  8 

Fig.  14  Optimum  &  values  corresponding  to  various  error  ratios  p, 
revised  formulation. 

discrepancy  measure  does  not  involve  the  measured  mode 
and  thus  cannot  be  used  by  itself.  In  fact,  it  is  necessary  that  the 
coefficient  amJ  that  minimizes  F\ „r)  also  provides  a  good 
representation  of  the  mode  shape  i/rln  .  This  requirement  will  be 


achieved  by  constraining  the  optimization  of  to  give  a  small 
representation  error;  that  is. 


l=i 


(30) 


To  impose  the  value  £„* ,  it  is  useful  to  first  determine  its  smallest 
value  Sm\  which  is  obtained  by  minimizing  E^'"J  [Eq.  (28)].  In  fact, 
Sm  is  the  representation  error  obtained  with  the  original  projection 
approach  to  determine  the  coefficients  amj.  Accordingly,  it  is 
proposed  here  to  set  £„’  =  p  £„  ,  where  p  is  a  user-selected  ratio 
greater  than  or  equal  to  1 . 

The  minimization  of  F m  [Eq.  (29)]  under  the  representation 
constraint  of  Eq.  (30)  can  be  achieved  by  the  Lagrange  multiplier 
(denoted  here  as  r)  approach,  i.e.,  with  the  unconstrained 
minimization  of 


;  (r)  _  ^  (a», 


■  +  T 


r„W 


I >™>./ 

1=1 


Jr) 


(3D 


Fig.  15  Maximum  value  of  the  log-likelihood  as  a  function  of  the  uncertainty  level  (A)  for  various  error  ratios  (p):  a)  p  =  1.31,  b)p  =  1.17,c)p  =  1.07, 
d)  p  =  1.03,  e)  p  =  1.016,  and  f)  p  =  1.008. 
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0  0.0022  -0.0152  0.0026  -0.0022  -0.0006 
-0.0022  0  -0.0342  0.0079  0.0018  0.0012 

0.0152  0.0342  0  -0.0688  0.0050  -0.0027 

-0.0026  -0.0079  0.0688  0  -0.0177  -0.0166 

0.0022  -0.0018  -0.0050  0.0177  0  -0.0640 

0.0006  -0.0012  0.0027  0.0166  0.0640  0 

0  0.0306  -0.0106  0.0035  -0.0002  0.0014  " 
-0.0306  0  -0.1002  0.0021  0.0008  -0.0003 

0.0106  0.1002  0  -0.0481  0.0010  0.0016 

-0.0035  -0.0021  0.0481  0  -0.0182  -0.0143 

0.0002  -0.0008  -0.0010  0.0182  0  -0.0319 

-0.0014  0.0003  -0.0016  0.0143  0.0319  0 

0  0.0344  -0.0126  0.0035  -0.0012  -0.0003' 
-0.0344  0  -0.1359  0.0006  0.0025  -0.0017 

0.0126  0.1359  0  -0.0655  0.0006  0.0016 

-0.0035  -0.0006  0.0655  0  -0.0410  -0.0207 

0.0012  -0.0025  -0.0006  0.0410  0  -0.0663 

0.0003  0.0017  -0.0016  0.0207  0.0663  0 


0  = 


0  = 


0  = 


0  -0.0025  -0.0052  0.0015  0.0006  0.0026 
0.0025  0  0.0926  0.0140  0.0022  0.0020 

0.0052  -0.0926  0  -0.0143  -0.0065  -0.0070 

-0.0015  -0.0140  0.0143  0  -0.0362  -0.0020 

-0.0006  -0.0022  0.0065  0.0362  0  -0.0230 

-0.0026  -0.0020  0.0070  0.0020  0.0230  0 

0  0.0072  -0.0081  0.0030  -0.0004  0.0003  ' 
-0.0072  0  -0.0189  0.0085  -0.0002  -0.0043 

0.0081  0.0189  0  -0.0565  0.0048  -0.0027 

-0.0030  -0.0085  0.0565  0  -0.0181  -0.0204 

0.0004  0.0002  -0.0048  0.0181  0  -0.0506 

-0.0003  0.0043  0.0027  0.0204  0.0506  0 

0  0.0276  -0.0086  0.0023  -0.0004  -0.0002 
-0.0276  0  -0.0778  0.0135  0.0071  -0.0021 

0.0086  0.0778  0  -0.0557  0.0026  -0.0025 

-0.0023  -0.0135  0.0557  0  -0.0072  -0.0036 

0.0004  -0.0071  -0.0026  0.0072  0  -0.0926 

0.0002  0.0021  0.0025  0.0036  0.0926  0 


Fig.  16  Values  of  obtained  in  the  six  tests.  Error  ratio  p  =  1.31. 


Differentiating  the  preceding  expression  with  respect  to  amj  leads 
to  the  system  of  equations 

[2  +  r$r$]  A«  =  t  +  £ Am  (32) 

where  2  is  the  diagonal  matrix  with  element  jj  equal  to  1  /<jA. 
Furthermore,  $  is  the  matrix  with  the  vector  0,  as  the  y'th  column,  and 
A{,n  and  A$  are  the  vectors  of  coefficients  <xmj,  and  amj. 

The  solution  of  the  linear  system  of  equations  of  Eq.  (32)  yields 
A{,n  in  terms  of  the  yet  unknown  Lagrange  multiplier  t.  Then, 
introducing  this  result  in  Eq.  (30)  leads  to  a  nonlinear  algebraic 
equation  for  r  in  terms  of  the  allowed  error  Once  solved,  the 
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Fig.  17  Values  of  the  mean  stiffness  matrix  K  obtained  with  the  error 
ratios  p  =  1.31. 


value  of  r  can  be  reintroduced  in  Am*  to  yield  the  final  expression  for 
the  coefficients  a,nj. 

It  should  be  noted  that  this  revised  approach  is  an  extension/ 
generalization  of  the  original  projection  method  described  before. 
Indeed,  if  one  sets  eln  =  £m\  the  previous  set  of  coefficients  is 
recovered  and  the  value  of  r  is  oc. 

Implicit  in  the  preceding  formulation  is  the  knowledge  of  the  mean 
matrix  K  and  the  uncertainty  level  X  that  is  necessary  for  the 
evaluation  of  amj  and  amj.  Certainly,  an  iterative  scheme  could  be 
developed  in  which  an  initial  estimate  of  K  and  X  is  used,  and  then  the 
process  could  be  updated  and  repeated  until  convergence.  In  this 
regard,  note  however  that  the  mean  values  amj  remain  very  close  to 
8mj  (the  Kronecker  symbol),  even  for  large  uncertainty  (small  X),  and 
thus  could  be  set  in  this  fashion.  Furthermore,  it  was  found  by 
experimenting  with  the  nonparametric  approach  that  the  standard 
deviations  amj  all  grow  approximately  proportionally  as  the  level  of 
uncertainty  increases,  except  for  the  diagonal  term  that  grows  much 
slower.  If  one  neglects  this  difference  of  the  diagonal  terms,  it  is 
found  that  the  solution  A„  is  independent  of  an  overall  scaling  of  the 
matrix  2,  which  would  only  be  reflected  in  the  value  of  the  Lagrange 
multiplier  r.  On  the  basis  of  this  discussion,  an  iterative  process  may 
not  be  necessary  and  amj  and  may  be  set  according  to  first 
estimates  of  the  mean  matrix  K  and  the  uncertainty  level  X. 

Notwithstanding  these  comments,  the  iterative  process  was 
implemented  here  to  fully  assess  the  methodology  (see  Fig.  13  for  a 
flowchart)  and  a  set  of  increasing  values  of  the  error  ratio  p  were 
considered.  As  indicated  in  Fig.  14,  the  allowance  of  a  slightly  larger 
error  permits  a  much  better  fit  by  the  nonparametric  methodology 
and  the  optimum  value  of  <5  that  decreases  rapidly  toward  the  value 
expected  from  a  fit  of  the  natural  frequencies  alone;  that  is,  as  p 
changes  from  1  to  1.31,  the  maximum  of  the  likelihood  function 
occurs  for  a  value  of  8,  which  decreases  from  0.6  to  0.05  (Fig.  15).  A 
study  of  the  resulting  values  (not  shown  here  for  brevity)  clearly 
demonstrates  that  the  outliers  highlighted  before  are  becoming 
smaller,  as  desired,  as  the  projection  error  is  allowed  to  increase 
slightly. 

These  results  demonstrate  that  the  revised  formulation  achieves 
exactly  the  stated  purpose,  and  thus  is  considered  here  as  the 
projection  approach  for  the  identification  of  the  maximum  entropy- 
based  uncertain  model.  But  how  should  the  error  ratio  be  set?  If  no 
specific  perspective  is  available  to  the  analyst,  the  error  should  be  set 
so  that  the  value  of  8  obtained  by  the  present  process  is  consistent 
with  the  one  that  can  be  estimated  from  the  natural  frequencies  alone, 
which  was  approximately  0.05  in  the  present  case.  Adopting  this 
criterion  leads  to  the  error  ratio  p  =  1.31,  for  which  the  final  values 
are  given  in  Fig.  16  and  lead  to  the  mean  stiffness  matrix  K  of 
Fig.  17. 

It  can  be  noted  on  the  two-dimensional  plots  of  Fig.  18  that  mode  6 
appears  to  have  a  larger  deviation  of  natural  frequencies  than  the 
other  ones.  On  the  basis  of  this  observation,  the  modified 
nonparametric  model  introduced  in  [3]  was  also  considered  with  a 
specific  variance  only  on  mode  6.  As  shown  in  Fig.  19,  the  ellipses 
corresponding  to  this  mode  have  increased  and  the  natural 
frequencies  of  mode  6  fit  within  these  ellipses  as  the  other  modes. 

All  results  described  previously  were  obtained  by  relying  on  the 
Gaussian  approximation  of  the  joint  probability  density  function  of 
the  random  vector  X  [see  Eq.  (9)]  to  build  the  likelihood  function  of 
Eq.  (6).  Nevertheless,  many  of  these  computations  were  also 
repeated  with  the  kernel  density  estimation  method  [1 1-13],  using 
the  Gaussian-based  optimum  solution  as  an  initial  condition,  and 
only  small  variations  in  the  parameters  (mean  stiffness  matrix  K  and 
uncertainty  parameter  X)  were  observed. 

VII.  Conclusions 

The  focus  of  this  investigation  was  on  the  formulation  and 
validation  of  a  methodology  for  the  estimation  of  an  uncertain  linear 
modal  model  of  a  structure  from  measurements  of  a  few  of  its  natural 
frequencies  and  mode  shapes  on  a  few  nominally  identical  samples 
of  the  structure.  The  basis  for  the  modal  model  is  composed  of  the 
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Fig.  19  Two-dimensional  joint  probability  density  functions  of  the  natural  frequencies  (freq)  and  values.  Error  ratio  p  =  1.31:  nonparametric 
model  of  [4]  for  the  last  mode. 
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modes  of  an  approximate  representation  of  the  structure,  e.g.,  a 
nonupdated  or  preliminary  FEM.  The  uncertainty  in  the  structure  was 
assumed  to  originate  from  stiffness  properties  with  the  stiffness 
matrix  represented  through  the  nonparametric  stochastic  modeling 
approach,  i.e.,  Eqs.  (1-3)  and  Fig.  1.  To  complete  the  modeling,  it  is 
thus  necessary  to  estimate  the  mean  stiffness  matrix  K  and  the 
parameter(s)  of  the  nonparametric  model,  which  was  achieved  using 
a  maximum  likelihood  strategy.  For  computational  expediency,  the 
joint  distribution  of  the  natural  frequencies  and  mode  shape  values 
necessary  in  this  strategy  was  approximated  by  a  joint  Gaussian 
probability  density  function. 

A  critical  step  in  this  effort  is  the  reduction  of  the  mode  shape  data 
to  infonnation  relevant  to  the  modal  model.  In  the  tests  considered 
here,  M  mode  shapes  of  every  sample  of  the  structure  would  be 
measured  at  a  large  number  nL  of  locations  that  likely  match  a  subset 
of  nodes  of  the  available  FEM.  However,  the  information  needed  for 
the  estimation  of  the  uncertain  modal  model  parameters  is  the 
representation  of  these  measured  modes  as  linear  combinations  of  the 
basis  functions.  This  task  was  first  accomplished  through  a  least- 
squares  projection  [see  Eq.  (5)]  of  the  measured  modes  on  the  basis. 
This  data  cannot  be  used  directly  for  the  estimation  of  the  uncertain 
model  parameters  because  of  the  lack  of  orthogonality  of  the 
measured  modes  with  respect  to  the  mass  matrix,  a  property  that  is 
automatically  satisfied  by  the  uncertain  modes  of  the  model.  The 
orthomormality  of  the  modes  was  seen  to  represent  constraints 
between  the  M  x  M  projection  coefficients  so  that  only  M  x  (M  — 
1  )/2  of  those  can  in  fact  be  considered  independent.  On  this  basis,  it 
was  proposed  to  reduce  further  the  modal  information  to  the  skew- 
symmetric  part  of  the  projection  coefficient  matrix;  see  Eqs.  (23-27). 

The  estimation  of  the  uncertain  model  parameters  from  such  data 
for  the  AFIT  joined  wing  led,  however,  to  an  uncertainty  level  much 
larger  than  could  be  expected  from  the  small  variability  of  the  natural 
frequencies.  A  detailed  analysis  of  this  issue  demonstrated  that  it 
arises  from  the  influence  of  noise  on  small  projection  coefficients. 
The  noise  effects  were  then  lessened  by  modifying  the  least-squares 
projection  of  the  measured  modes  on  the  basis:  allowing  a  slightly 
larger  error  in  the  projection  but  constraining  a  better  fit  with  the 
uncertain  modal  model  structure;  see  Eqs.  (29-32)  and  Fig.  13.  The 
estimation  of  the  uncertain  modal  model  parameters,  i.e.,  the  mean 
stiffness  matrix  K  and  the  parameter(s)  of  the  nonparametric 
representation,  was  then  achieved  successfully;  see  Figs.  17-19.  The 
methodology  proposed  here  is  applicable  to  a  broad  class  of  linear 
structures  and  data,  e.g.,  accelerometer  vs  laser  vibrometer,  dense  vs 
sparse  measurements,  etc. 
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